Reaction-diffusion front crossing a local defect 
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The interaction of a Zeldovich reaction-diffusion front with a localized defect is studied numerically 
and analytically. For the analysis, we start from conservation laws and develop simple collective 
variable ordinary differential equations for the front position and width. Their solutions are in good 
agreement with the solutions of the full problem. Finally using this reduced model, we explain the 
pinning of the front on a large defect and obtain a quantitative criterion. 

PACS numbers: 



I. INTRODUCTION 

Since the pioneering work of Zeldovich et al [l| re- 
action diffusion equations have been important mod- 
els to describe the combustion of solid materials. The 
so-called combustion front connects two equilibria, the 
unburnt solid on one side and the burnt solid on the 
other. These models also describe in some limit the 
propagation of a nerve impulse in a neuron. See 0] for 
a review of these applications. FisherQ in 1937 with 
his model of gene propagation opened a new area of 
application. In all these systems there is one compo- 
nent, a concentration which diffuses and a reaction term, 
which is cubic for the Zeldovich model and quadratic 
for the Fisher equation. From another point of view, 
these one variable systems can be seen as a reduction 
of more general models like the Susceptible-Infected- 
Recovered (SIR) model of Kermac-McKendrick [|| de- 
scribing the propagation of epidemics. An easy way to see 
this is to start from the Susceptible-Infected-Recovered 
(SIR) three-component quadratic model of Kermack- 
McKendrick Q. We assume R = so that S = N - I 
where N is the total number. The quadratic terms are 
then S I = (N — 1)1. Therefore the Fisher or Zeldovich 
equations can be considered as the simplest model of the 
propagation of an epidemic. 

In many cases the reaction term is non-uniform in 
space. The powder in a combustion tube can present 
some defects along the tube. A sudden enlargement of 
an axon is known to block nerve impulse as discussed in 
Q. This effect was analyzed by Chapuisat et al @. Ge- 
ographical factors like rivers or mountains slow down an 
epidemic front. This can be seen on the records of the 
plague epidemic in the 13th century. See for example [6] 
and Q on this example. These factors also need to be 
considered when modeling the diffusion of certain genes 



[8| . Some defects are also time dependant like seasonality 
effects. In a number of situations, the defect can be con- 
sidered as localized in a given region of space. This is the 
case for a river in the example of the propagation of an 
epidemic given above. Such a local defect will enable to 
reduce the dimensionality of the problem. For example 
a line defect in 2D will separate the plane in two so that 
one can consider only the propagation along the ID line 
normal to the defect line. In the following article on the 
Zeldovich model we study the interaction of a ID front 
with a localized defect. We consider a ID model because 
it is simple, one can do analysis and there is an exact 
front solution if the reaction term is a third degree poly- 
nomial. As mentioned above even if the problem is 2D 
or 3D, a localized defect reduces the geometry to ID: be- 
fore the defect and after the defect. This approximation 
makes sense of course for a narrow tube. For this model, 
we show that the front can be stopped by a large enough 
defect. We introduce an approximate analysis based on 
conservation laws which gives ordinary differential equa- 
tions for the front position and width. The solutions of 
these simple models are in good agreement with the solu- 
tion of the partial differential equation. These collective 
variable equations easily lead to the solution of the in- 
verse problem of determining the defect from the motion 
of the front. 

The article is organized as follows. Section HT1 presents a 
preliminary analysis of the model. Numerical solutions 
are shown and analyzed in section IIIII Section IIVI intro- 
duces the collective variable equations whose solutions 
are compared to the solutions of the full equation in sec- 
tion El We conclude in section IVlTl 



II. PRELIMINARY ANALYSIS 



A general reaction-diffusion equation in an heteroge- 
neous medium is 
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U t = U xx + s(x)R(u) 



(1) 
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We consider that R is a third degree polynomial with 
three real roots so that the model ([!]) is written 



iii 



+ s(x)u(l — u)(u — a) 



(2) 



with < a < 1. The homogeneous states are u* = 0, 1, a. 
Only the first two are stable [2[ so we will consider fronts 
connecting to 1. If s is homogeneous, the front solution 
can be calculated by assuming that it is a traveling wave 
u{z) = u{x — ct) and assuming that du/dz = Ku(l — u) 
where K is a, constant @, 

One obtains the kink exact solution 



u(x. t) 



1 



l + exp(± v /f(a;-ct)) 



where the speed is 



C = ±J-(l-2a). 



(3) 



(4) 



Depending on the ± sign we have a increasing (resp. de- 
creasing) kink going from (resp. 1) as x — > — oo to 1 
(resp. 0) as x — > +oo for the + (resp. — ) sign. The 
width of the kink is 



(5) 



It is inversely proportional to the inhomogeneity so a 
large s corresponds to a very sharp and fast front. 



III. NUMERICAL ANALYSIS FOR A 
LOCALIZED DEFECT s(x) 



We will consider two main localized defects : a gaus- 
sian defect and a tanh defect. There are two length scales 
in the problem, the kink width w given by ([5]) in the ho- 
mogeneous case and the characteristic length d of the de- 
fect. For a given front, the defect will be wide or narrow 
depending on the ratio of these two length scales. 



A. Gaussian defect 



Specifically the first class of defect is given by 



s(x) = sq + Si exp 



2d 



(7) 



When the defect varies on a length scale large compared 
with the width of the kink, the front moves adiabatically. 
An example is shown in figure [TJ The left panel shows 
the initial front together with the gaussian defect as a 
function of position x. The right panel shows the speed 
x' and width w of the kink as a function of the kink 
position xo- The defect s(xo) is also reported in this 
graph. The speed is computed using a centered difference 
approximation from the time-series xo(t). One sees that 
the kink accelerates and gets steeper as it runs into the 
defect. 



When s(x) varies there is only a local balance be- 
tween the reaction term R and the diffusion term u xx . 
The equation does not have an explicit solution so we 
study it numerically. Specifically we use the method of 
lines where we discretize space by the finite difference 
method. The time advance is given by an ordinary differ- 
ential equation solver. All the details are given in the ap- 
pendixjA] For simplicity we chose the parameter a = 0.3 
throughout the article. 

The fronts are stable for many different values of the 
localized defect s(x). Precisely, the fronts continue to 
exist with a speed and a width that change. Therefore it 
is reasonable to fit the solution at each time using a least 
square procedure by a kink 



u k (x,t) = 



1 



1 + exp 



/ x-xp(t) \ 

\ Mt) ) 



(6) 



where the time dependence is only through the collective 
coordinates, the kink position xq and width w. The de- 
tails of the fitting procedure are given in the appendix [A"l 
The important fact is that the normalized least-square er- 
ror is small (< 10~ 4 ) in all the numerical results. There- 
fore the fit is very good. 
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FIG. 1: Adiabatic motion of a front. The left panel shows 
the initial front u(x,t = 0) together with the defect s(x). 
The right panel shows the speed x'q and width w of the kink 
as a function of the kink position xq. The parameters are 
so = 0, 3, si = 0, 6 and d = 30. 



Now let us consider the other limiting case, i.e. when 
the defect is narrower than the kink as shown in the left 
panel of figure [2j Here again the kink is able to breach the 
obstacle although it slows down before it and accelerates 
past it as shown in the right panel of figure [5] Notice 
also the modulation of the width. 



FIG. 2: Motion of a front as it hits a narrow defect. The 
left panel shows the initial front u(x, t — 0) together with the 
defect s(x). The right panel shows the speed x'o and width 
w of the kink as a function of the kink position xq. The 
parameters are so = 0, 3, si = 0, 6 and d = 0, 3. 

If the amplitude of the defect is larger, the kink can 
be stopped as shown in figure [3] The left panel presents 
a defect with an amplitude s\ = 7. The right panel 
shows the width and the speed of the front. The latter 
goes to zero and the width remains stationary. To show 
that this effect is intrinsic we varied systematically the 
spatial resolution of the computation from N = 800 to 
N = 12000. For all these runs the front stopped at the 
same position xq ~ —2, 73, with a width w ~ 1, 80. This 
phenomenon is also called "pinning of the front" 0. We 
will analyze it in detail below. 
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FIG. 4: Interaction of the front with a wide "tanh" defect. 
The parameters are s; = 0, 3, s r = 1 (jump of 0, 7) and d = 10. 

When the defect is narrow, the speed and width be- 
come non monotonic. The results are shown in figure [5] 
Again the results have been checked by varying system- 
atically the spatial resolution. 
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FIG. 5: Interaction of the front with a narrow "tanh" defect. 
The parameters are s; = 0, 3, s r = 1 (jump of 0, 7) and d = 
0,1. 

As previously the dip in the velocity indicates that 
increasing the amplitude we will obtain the pinning of 
the front. This indeed happens. We will show in section 
IVII that we can infer a simple criterion for the pinning of 
the front. 



IV. COLLECTIVE COORDINATE ANALYSIS 



FIG. 3: Pinning of the front when it hits a large defect. The 
parameters are so = 0, 3, si = 7 and d = 0, 3. 



B. Tanh defect 



The numerical results in the previous section can be 
understood in a simple way. Assume that the solution 
is a function of the reduced variable z — x — ct where 
the speed c is a slowly varying function of time. Then 
reporting u into the partial differential equation (JXJ) and 
integrating we get 



We now consider a different type of defect which con- 
nects two different values s; on the left and s r on the 
right. 

S (x) = S; + ^ y ^(l + tanh^) (8) 

We observe similar results as for the gaussian defect. For 
a "wide" defect the front moves adiabatically as shown 
in figure HI It speeds up and its width decreases as s 
increases from left to right. 



/+oo 
U'(x — ct)dx 
-co 

/+oo /- + oo 

U"(x-ct)dx+ / s(x)R(u(x,t))dx (9) 
-oo J — oo 

Consider a front going from 1 (to the left) to 0. Since 
the front is flat away from the defect we have 

J U"(z)dz = 0, J U'(z)dz = -1, 



4 



so that we obtain 



+00 



s(x)R(u(x, t))dx. 



(10) 



This result shows that the speed decreases (resp. in- 
creases) when the front reaches a region where s is smaller 
(resp. larger). This is not completely correct because the 
width of the front also depends on s as shown in ((3J). 

A more general approach is then to assume that the 
front keeps its general form but that its center xo and 
width w become modulated 



u{x, t) = U 



x {t) 



w(t) 



U(z) 



(11) 



At this point we do not specify the form of U(z). This 
method is often called "variational approach" in the 
context of Hamiltonian partial differential equations Q 
which are derived from a Lagrangian density. Here how- 
ever there is no variational principle so we need to gen- 
erate two conservation laws to obtain the evolution of 
xq and w. The first one is the original partial differen- 
tial equation (p}. The second conservation law can be 
obtained by multiplying the equation by u 



uu t — uu xx + s(x)uR(u). 



(12) 



In principle we could have also multiplied by u x but this 
would give the wrong evolution if s is singular. After 
replacing the expression (|11|) into the two partial differ- 
ential equations ([T|). (IT2l . we get the evolutions of the 
front position xq and width w 



J J +w J s ( wz + x o)R = 0, 

x' J UU'+w' J UU'z+w J s{wz+x )UR = i J U' 2 , 

(13) 

where the terms fU', JU'z,fUU', JlIU'z are inte- 
grals with respect to the variable z and therefore are 
numbers. The integrals 



s{wz + Xq)R 



s(wz + xq)UR = 



s(wz + xo)R(U(z))dz, 



+00 



s(wz + xo)U(z)R(U(z))dz, 



depend on xq and w. They yield the source terms for 
the differential equations. All the details are given in 
Appendix B. These ordinary differential equations are 
very general and can be transposed to different types of 
nonlinearities. The only assumptions are that the defect 
s(x) is localized and that the front keeps its functional 
profile. 



We consider the simplest ansatz jBJ which is suitable 
for the Zeldovich equation because it is an exact solution 
when s is constant. Then ((15)) reduce to 



+00 



s(wz + xo)R(U(z))dz, 



1 

3w 



+00 



s(wz + x )(l - 2U(z))R(U(z))dz. 

Note that we still have made no assumptions on the non- 
linearity R. 

We will consider three main situations, a defect that 
varies on a scale longer than the width of the front and 
two sharp defects represented respectively as a Dirac dis- 
tribution and a Heaviside function. For a defect that 
varies on a scale longer than w, s(xq + wz) w s(xq) so 
that s goes out of the integral, modifying the equations 
subsequently 



1 - 2a 
2 

3u> 6 



ws(xo), 



1 w 
w = -ts(xo). 



(15) 



These equations justify the notion of a local speed and 
width of the front. When s(x) = s is constant, we re- 
cover the results ((4|, ((5]). Another remark is that here 
the defect can be extended. The only condition that we 
require is that the front has reached it's equilibrium state 
before reaching the defect region. 

Note that the approximation can be sharpened by use 
of a Taylor expansion : s(xq + wz) w s(xq) + wzs'(xq). 
Then the equations take the form 



. 1 w . . (1 — 2a)w 2 ,. 
W = ^-6 3{Xo) + 2 S( * 0) - 



(16) 



We now assume a sharp "bump-like" defect described 
by s(x) — a + f3S(x) where 8(x) is the Dirac distribu- 
tion centered at x = 0. The system of equations is now 
reduced to 



in = aw 



1 - 2a 



+ /3RIU 



-Xq 

w 



The source terms are given by 



R{U{z)) = -e 



(17) 



— 1 + a + ae z 



(l-2U(z))R(U(z)) = e/(l-e*) 



They are shown in Fig. |B1 



(l + e z ) 3 ' 

1 + a + ae z 



(18) 
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FIG. 6: Plot of the source terms R(U(— z)) (left panel) and 
(1 — 2U{— z)) R(U( — z)) (right panel) for a Dirac distribution 
defect as a function of the reduced variable z, for a chosen 
parameter a = 0.3. 

Note how the two sources terms R(U(—z)) and 
(1 — 2U{— z)) R(U(—z)) are not symmetric. Because 
R(U(—z)) is largely positive a defect will cause an ac- 
celeration of the kink for /3 > 0. But because there is 
also a negative part, the pinning of the front becomes 
possible both for (3 > and (3 < 0, see section PVTI for de- 
tails. Notice also how (1 — 2Z7(— zj) R(U(—z)) is largely 
negative so that the width of the front decreases as it hits 
the defect for /? > 0. 

The other sharp defect that we study is s(x) = a + 
j3H{x) where H (x) is the Heaviside function. The system 
of equations is now reduced to 



The general remarks on the acceleration and width re- 
duction of the front still apply. 

To summarize, we have obtained fairly simple ordinary 
differential equations for the evolution of the front posi- 
tion xq and width w for the Zeldovich equation. We will 
see in the next section that these equations yield very 
good approximations to the solution of the partial differ- 
ential equation. 



V. COMPARISON BETWEEN THE FULL 
MODEL AND THE REDUCED MODEL 

To establish the validity of the reduced model, it is 
important to compare its solutions to the ones of the 
partial differential equation. As discussed in the previous 
section, we classify the defects s(x) as wide or narrow 
depending whether w/d <C 1 or w/d> 1 where w is the 
initial width of the front and d is the width of the defect 
as defined in formulas ([7]) and ©. 



w 



aw 
1 

3w 



1 - 2a 



(3w 



W 



r+oo 
Pw I (1 



Here the source terms arc 



R(U(z))dz, 

- 2U{z))R{U{z))dz. 

(19) 



A. Adiabatic case 

When the defect is wide, equations (fTSj) give evolutions 
of the front position xq and width w that are close to the 
fits obtained from the solutions of the partial differential 
equation . The plots are shown in Fig. [5] for the gaussian 
defect of Fig. Q]and the tanh defect of Fig. @] 



+oo 



(1 



R(U(z))dz 
2U{z)) R{U{z))dz 



They are plotted in Fig. [71 






FIG. 8: Plots (xo, w) for the partial differential equation ([2]) in 
continuous line (red online) and the simple collective variable 
equations (115[l in dashed line (blue online). The left panel 
(resp. right panel) corresponds to a gaussian (resp. tanh) 
defect. The defects are shown at the bottom of each panel. 



FIG. 7: Plot of the source terms f+™ R(U(z))dz (left panel) 

and f_ °° (1 — 2U(z)) R(U ( z))dz (right panel) for a Heaviside 
distribution defect as a function of the reduced variable y, for 
a chosen parameter a — 0.3. 



The collective variable estimates can be improved by 
calculating numerically the integrals in equation (|13[) . 
We compare in Fig. |H] the results for the partial differ- 
ential equation (f5]) and for the solutions of (IT51) for the 
"tanh" defect. As can be seen the agreement is excellent. 
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FIG. 9: Plot (xo, w) for the partial differential equation ([2]) in 
continuous line (red online) and the collective variable equa- 
tions (|13l) in dashed line (blue online) for the "tanh" defect. 



B. Narrow defect: pinning 

We now consider that the defect width d is smaller than 
the front width w. Then the collective variable equations 
can be reduced as shown in the previous section. We 
approximate the gaussian s(x) = sq + si exp f by a 
Dirac delta function s(x) — a + (3S(x). We choose a = sq 
and (3 — siV^nd so that 



^ = jL SieXP {-2d) dx 



The integrals associated to the defect are then equal. We 
tested the validity of this approximation and found that 
the solutions agree to about 1 % when w > lOd. 

For such narrow defects, the collective variable equa- 
tions (fl7|) are less accurate than for a wide defect. They 
do provide however the qualitative behavior, in partic- 
ular the pinning of the front. Fig. [10] shows the plots 
for a gaussian defect drawn at the bottom. The value 
of s at infinity is sq = 0,3. The left panel corresponds 
to a small amplitude s\ = 0, 6. The right panel is for 
a much larger amplitude si = 5 causing the pinning of 
the front. For this particular plot we show the speed x' 
computed using a centered difference for both the partial 
differential equation and the collective variables. They 
are multiplied by 10 in the plot for clarity. The defect 
has been divided by 10. 



FIG. 10: Plots (xo,w) for the partial differential equation 
([2]) in continuous line (red online) and the collective variable 
equations (|17[) in dashed line (blue online) for a narrow gaus- 
sian defect shown at the bottom of the plots. The left panel 
corresponds to si = 0, 6 and d — 0,3. The right panel corre- 
sponds to si = 5. There the speed x' is also reported. 



The partial differential equation and the collective vari- 
ables agree well for the width w as a function of xo . The 
speed x' is not so well approximated but it goes to zero 
for a pinning position xq. 

For narrow "tanh" defects, we reduce the collective 
variable equations to the ones for a Heaviside defect (fT9]) 
as shown above. The agreement between the curves 
(xq,w) for the partial differential equation @ and the 
collective variable equations (fT§|) is good as shown in Fig. 
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FIG. 11: Plot (xo, w) for the partial differential equation ([2]) in 
continuous line (red online) and the collective variable equa- 
tions (|19[) in dashed line (blue online) for a narrow "tanh" 
defect shown at the bottom of the plot. The parameters are 
si — 0, 3, s r = 1 -corresponding to a jump of 0, 7- and d — 0, 1. 



As we have seen in section [TTT1 the front can get pinned 
when si is large enough and this is predicted also by 
the the collective variable model. Such an example is 
shown in Fig. [12] for a large and narrow defect where 
s r = 8, d = 0, 1 shown scaled by 0.1 at the bottom of the 
plot. As previously the speed x' estimated using finite 
differences has been plotted in the graph to show pinning. 
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FIG. 12: Plot (xo, w) for the partial differential equation ([2]) in 
continuous line (red online) and the collective variable equa- 
tions (|19[) in dashed line (blue online) for a narrow "tanh" 
defect shown at the bottom of the plot. The speed x'q is also 
reported using the same color scheme. 



VI. ESTIMATES FOR FRONT PINNING AND 
DEFECT TOPOGRAPHY 



We now illustrate how the reduced model (fT4"f can be 
used. Its main advantage is that the parameters ap- 
pear explicitly and therefore their influence can be un- 
derstood. A first direct application is a simple criterion 
for the front pinning. 

From the equation for the Dirac delta function defect, 
we can obtain a rough estimate of the strength /3 of the 
defect necessary to stop the front. The evolution of the 
front from (fT7|) shows that the front can stop x' — 0, w' — 
when 
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1 - 2U 



-x 



-x 







R U 



-x 



(21) 

Therefore the front can stop when the two terms on the 
right hand side of the first equation balance each other. 
In other words we need that 



aw 



2a 



^ -/3min(R). 



So there exist a threshold for the pinning : 
1 - 2a" 



= — /3 C min(i?). 



(22) 



The minimum of R(U(z)) as a function of z and its ar- 
gument z m i n can be computed exactly, 



= In 



1 



i(R) = R(U(z mm )) = 



a 2 (l + r)(a + r) 
(1 + a + r) 3 ' 



(23) 



(24) 



where r = yfl — a + a 2 . Plugging f3 = (3 C into the second 
equation of (|2"Tj) and assuming — z m i n , we get an 
expression for the width at pinning w. Using this we also 
get the front stopping position a;o- For the critical /3, the 
position and width at pinning are given by 



Xo 



l + 3(l-2a)(i=2±£ 
2 1 



l + 3(l-2a)fe£ 



(25) 



The critical /3 is obtained from (|22)) with the above value 
of w. It reads 



(1 - 2o)(l + a + r f 



a 2(l + r )( a + r )Jl + 3(l-2a)(i^ 



These values are in good agreement with the values ob- 
tained numerically in the previous section, see Fig. |3] 
Precisely, with a = 0, 3, and a — 0,3 too, we find /3 C ~ 
5,88, and the parameters at pinning xq ps —3,47, w ~ 
1, 89. Numerically, we found j3 c ps 6, 18, xo ~ —2, 73, w pa 
1,80. This (3 C corresponds to Si pa 7,80 because we set 
the width of the defect d — 0, 1. Such a criterion can also 
be infered in the case of a Heaviside defect. A pinning of 
the front is also possible with f3 < (as long as (3 > — a 
so that s(x) remains positive). 

Another direct application of the collective variable 
model is that we can obtain the defect topography s(x) 
from the observation of the front position xo and width 
w. We illustrate this on the example of a wide gaus- 
sian defect of the form jJJ with (|2|) and the parameters 
so = 0, 6, s± = 0, 3, d = 10. The position and width of 
the front are estimated using the least square fit on the 
solution of the partial differential equation ([2]) . From the 
collective variable equations in the adiabatic case (jT5j) we 
get 



s(x ) 



'o 



1 — 2a w 



(27) 



Using a centered difference approximation for the time 
derivative, we obtain an estimate of s(xo)- This estimate 
is compared to the "real" s(x) in Fig. Q21 
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FIG. 13: Defect topography s(x) estimation from the solution 
of the partial differential equation. The estimated defect to- 
pography s(x) using (|27|l is shown in dashed line (blue online) 
while the "real" s(x) is shown in continuous line (red online). 

As can be seen the agreement is very good. This ap- 
proach can then be extended to other types of defects for 
physical or biological applications. 

These two examples illustrate the power of these re- 
duced models. The role of the parameters is very easy 



to understand. One can easily solve the inverse problem 
of estimating these parameters from measurements. An- 
other extension of this could be to control the front using 
these reduced equations. 



VII. CONCLUSION 

We have analyzed numerically the interaction of Zel- 
dovich reaction-diffusion front with a localized defect. 
The stability of the front for different types of defects 
suggested that it has the form of a generalized traveling 
wave with a time dependant position and width. Us- 
ing conservation laws we obtained ordinary differential 
equations for these collective variables. We further re- 
duced these models for the cases of an adiabatic defector 
a sharp "gaussian" or "tanh" defect. For these three 
cases the position and width obtained by fitting the nu- 
merical solution agree very well with the solutions of the 
collective variable equations. Finally we illustrated how 
these reduced models can be used to predict the pinning 
of the front on a large defect or to estimate the defect to- 
pography from a time-series of front positions and width. 
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Appendix A: Numerical procedures 

The basis of the method is to discretize the spatial 
part of the operator and keep the temporal part as such. 
We thereby transform the partial differential equation 
into a system of ordinary differential equations. This 
method allows to increase the precision of the approxi- 
mation in time and space independently and easily. We 
choose as space discretisation the finite volume approx- 
imation where the operator is integrated over reference 
volumes. The value of the function is assumed constant 
in each volume. As solver for the system of differential 
equations, we use the Runge-Kutta method of order 4- 
5 introduced by Dormand and Prince which enables to 
control the local error by varying the time-step. This 
has been implemented as the Fortran code DOPRI5 by 
Hairer and Norsett [Icj . 

We introduce reference volumes V k (interval) whose 
centers are at discretisation points x k = x m in + h/2 + 
(k - l)h, 1 < k ^ n with h = (x max - x min )/n 



X k + X k+ l + x k 



, 1 < k < n. 



For a fixed t, we assume u(x,t) to be constant on each 
volume Vk, u(xk,t) — u k . Integrating over V k , we obtain 
for 1 < k < n 



Uk 



h 2 



s k R(u k ) 



(Al) 



The boundary points k = 1, n are such that one satisfies 



9 



the homogeneous Neumann boundary conditions at x 

Xjnin i Xjnax • ^Ve have 



U2 ~ Ui 

ui = — p R(ui)si, 



-U n + U n -i 
Un = ~ h R(U n )S n . 



h 2 



(A2) 



(A3) 



The fit of the solutions using the front profile is done 
using a least square method. For each time t we introduce 
an error function 



1 " 

E(x , w) = - y^(ui(t) - u k {x l , t)Y 



(A4) 



»=i 



where the ttj are the values of u calculated numerically 
and 



u k {x,t) 



1 



1 + exp 



x-x (t) 
w(t) 



is the exact kink solution. The expression E is mini- 
mized usin g th e Polak-Ribiere combination of line min- 
imisations [llj . The fit of the solutions is done on the 
region such that 0.01 < u k < 0.99. The initial guesses 
for xq and w are estimated from u s» 0.25 and u ps 0.75. 
The method converges in about 20 iterations and the 
value of the energy is small, min E sC 10~ 4 . Therefore 
the fit is good. 



get 



x' . ( x - x \ ^ f + °° w'(x - s ) TTl fx- x 



>-U' 
to V iv 

r+oo 



dx- 



-U' 



8 ( X )R I u ( ) ) dx (B5) 



We compute the integrals by making the change of vari- 
ables x = w z + Xq 



/+oo r + oo 

U'(z)dz — w' I U'(z)zdz 
-oo J — oo 

/+oo 
s(wz + xo)R(U(z))dz 
-oo 

Assuming that the front goes from 1 to and that U'(z) 
is even we get the final result 

/+oo r+oo 
U'(z)dz + to' I U'(z)zdz 
-oo J —oo 

r+oo 

+ w s{wz + x )R(U(z))dz = (B6) 



The second conservation law (1121) will yield the evolu- 
tion of w. Proceeding as above we get 



Appendix B: Derivation of the collective variable 
equations 

The traveling wave is 

u(x,t) = ui^-^j^U(z), (Bl) 

so that we have the relations for the partial derivatives 
If, w'(x-x )\ ,fx-x ' 



Ui 



1 , f x - s 



If 



to V to 



Hxx r> U 



1 ttII f X- X 



, (B2) 
(B3) 
(B4) 



The first equation is obtained by integrating (JT]) with 
respect to s 



/+oo r-too 
u t dx = Kj+^-oo + / s(x)R(u)dx 
-oo J —oo 

The term [iix]j^-co i s zero because the front is flat away 
from the defect. Introducing the partial derivatives we 



/+oo r+oo 
U(z)U'(z)dz + w' U(z)U'(z)zdz 
-oo J —oo 

+ oo i r+oo 

s(wz+x a )U(z)R(U(z))dz = - U' 2 {z)dz 

-oo J —oo 

(B7) 

Note that we only assumed U — 1 at — oo and U = 
at oo and that U' is even. These assumptions are very 
general. In particular we have made no restrictions on 
the reaction term R. 

For the Zeldovich reaction term R(u) = u(l — u)(u — a) 
it is natural to assume that 



U{z) = 



1 



1 + exp z 



Then the integrals not involving s can be computed and 
we obtain the final result. 



+oo 



o = to / s(wz + xo)R(U(z))dz 
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The interaction of a Zeldovich reaction-diffusion front with a localized defect is studied numerically 
and analytically. For the analysis, we start from conservation laws and develop simple collective 
variable ordinary differential equations for the front position and width. Their solutions are in good 
agreement with the solutions of the full problem. Finally using this reduced model, we explain the 
pinning of the front on a large defect and obtain a quantitative criterion. 

PACS numbers: 



I. INTRODUCTION 

Since the pioneering work of Zeldovich et al [l| re- 
action diffusion equations have been important mod- 
els to describe the combustion of solid materials. The 
so-called combustion front connects two equilibria, the 
unburnt solid on one side and the burnt solid on the 
other. These models also describe in some limit the 
propagation of a nerve impulse in a neuron. See 0] for 
a review of these applications. FisherQ in 1937 with 
his model of gene propagation opened a new area of 
application. In all these systems there is one compo- 
nent, a concentration which diffuses and a reaction term, 
which is cubic for the Zeldovich model and quadratic 
for the Fisher equation. From another point of view, 
these one variable systems can be seen as a reduction 
of more general models like the Susceptible-Infected- 
Recovered (SIR) model of Kermac-McKendrick [|| de- 
scribing the propagation of epidemics. An easy way to see 
this is to start from the Susceptible-Infected-Recovered 
(SIR) three-component quadratic model of Kermack- 
McKendrick Q. We assume R = so that S = N - I 
where N is the total number. The quadratic terms are 
then S I = (N — 1)1. Therefore the Fisher or Zeldovich 
equations can be considered as the simplest model of the 
propagation of an epidemic. 

In many cases the reaction term is non-uniform in 
space. The powder in a combustion tube can present 
some defects along the tube. A sudden enlargment of an 
axon is known to block nerve impulse as discussed in Q . 
This effect was analyzed by Chapuisat et al Q. Geo- 
graphical factors like rivers or mountains slow down an 
epidemic front. This can be seen on the records of the 
plague epidemic in the 13th century. See for example [6] 
and Q on this example. These factors also need to be 
considered when modeling the diffusion of certain genes 



[8| . Some defects are also time dependant like seasonality 
effects. In a number of situations, the defect can be con- 
sidered as localized in a given region of space. This is the 
case for a river in the example of the propagation of an 
epidemic given above. Such a local defect will enable to 
reduce the dimensionality of the problem. For example 
a line defect in 2D will separate the plane in two so that 
one can consider only the propagation along the ID line 
normal to the defect line. In the following article on the 
Zeldovich model we study the interaction of a ID front 
with a localized defect. We consider a ID model because 
it is simple, one can do analysis and there is an exact 
front solution if the reaction term is a third degree poly- 
nomial. As mentioned above even if the problem is 2D 
or 3D, a localized defect reduces the geometry to ID: be- 
fore the defect and after the defect. This approximation 
makes sense of course for a narrow tube. For this model, 
we show that the front can be stopped by a large enough 
defect. We introduce an approximate analysis based on 
conservation laws which gives ordinary differential equa- 
tions for the front position and width. The solutions of 
these simple models are in good agreement with the solu- 
tion of the partial differential equation. These collective 
variable equations easily lead to the solution of the in- 
verse problem of determining the defect from the motion 
of the front. 

The article is organized as follows. Section HT1 presents a 
preliminary analysis of the model. Numerical solutions 
are shown and analyzed in section IIIII Section IIVI intro- 
duces the collective variable equations whose solutions 
are compared to the solutions of the full equation in sec- 
tion El We conclude in section IVlTl 



II. PRELIMINARY ANALYSIS 



A general reaction-diffusion equation in an heteroge- 
neous medium is 
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U t = U xx + s(x)R(u) 



(1) 



2 



We consider that R is a third degree polynomial with 
three real roots so that the model ([!]) is written 



iii 



+ s(x)u(l — u)(u — a) 



(2) 



with < a < 1. The homogeneous states are u* = 0, 1, a. 
Only the first two are stable [2[ so we will consider fronts 
connecting to 1. If s is homogeneous, the front solution 
can be calculated by assuming that it is a traveling wave 
u{z) = u{x — ct) and assuming that du/dz = Ku(l — u) 
where K is a, constant @, 

One obtains the kink exact solution 



u(x. t) 



1 



l + exp(± v /f(a;-ct)) 



where the speed is 



C = ±J-(l-2a). 



(3) 



(4) 



Depending on the ± sign we have a increasing (resp. de- 
creasing) kink going from (resp. 1) as x — > — oo to 1 
(resp. 0) as x — > +oo for the + (resp. — ) sign. The 
width of the kink is 



(5) 



It is inversely proportional to the inhomogeneity so a 
large s corresponds to a very sharp and fast front. 



III. NUMERICAL ANALYSIS FOR A 
LOCALIZED DEFECT s(x) 



We will consider two main localized defects : a gaus- 
sian defect and a tanh defect. There are two length scales 
in the problem, the kink width w given by ([5]) in the ho- 
mogeneous case and the characteristic length d of the de- 
fect. For a given front, the defect will be wide or narrow 
depending on the ratio of these two length scales. 



A. Gaussian defect 



Specifically the first class of defect is given by 



s(x) = sq + Si exp 



2d 



(7) 



When the defect varies on a length scale large compared 
with the width of the kink, the front moves adiabatically. 
An example is shown in figure [TJ The left panel shows 
the initial front together with the gaussian defect as a 
function of position x. The right panel shows the speed 
x' and width w of the kink as a function of the kink 
position xo- The defect s(xo) is also reported in this 
graph. The speed is computed using a centered difference 
approximation from the time-series xo(t). One sees that 
the kink accelerates and gets steeper as it runs into the 
defect. 



When s(x) varies there is only a local balance be- 
tween the reaction term R and the diffusion term u xx . 
The equation does not have an explicit solution so we 
study it numerically. Specifically we use the method of 
lines where we discretize space by the finite difference 
method. The time advance is given by an ordinary differ- 
ential equation solver. All the details are given in the ap- 
pendixjA] For simplicity we chose the parameter a = 0.3 
throughout the article. 

The fronts are stable for many different values of the 
localized defect s(x). Precisely, the fronts continue to 
exist with a speed and a width that change. Therefore it 
is reasonable to fit the solution at each time using a least 
square procedure by a kink 



u k (x,t) = 



1 



1 + exp 



/ X~Xg(t) \ 

\ Mt) ) 



(6) 



where the time dependance is only through the collective 
coordinates, the kink position xq and width w. The de- 
tails of the fitting procedure are given in the appendix [A"l 
The important fact is that the normalized least-square er- 
ror is small (< 10~ 4 ) in all the numerical results. There- 
fore the fit is very good. 
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FIG. 1: Adiabatic motion of a front. The left panel shows 
the initial front u(x,t = 0) together with the defect s(x). 
The rigth panel shows the speed x' and width w of the kink 
as a function of the kink position xq. The parameters are 
so = 0, 3, si = 0, 6 and d = 30. 



Now let us consider the other limiting case, i.e. when 
the defect is narrower than the kink as shown in the left 
panel of figure [2j Here again the kink is able to breach the 
obstacle although it slows down before it and accelerates 
past it as shown in the right panel of figure [5] Notice 
also the modulation of the width. 
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FIG. 2: Motion of a front as it hits a narrow defect. The 
left panel shows the initial front u(x, t — 0) together with the 
defect s(x). The right panel shows the speed x'o and width 
w of the kink as a function of the kink position xq. The 
parameters are so = 0, 3, si = 0, 6 and d = 0, 3. 



If the amplitude of the defect is larger, the kink can 
be stopped as shown in figure [3] The left panel presents 
a defect with an amplitude s\ = 7. The right panel 
shows the width and the speed of the front. The latter 
goes to zero and the width remains stationary. To show 
that this effect is intrinsic we varied systematically the 
spatial resolution of the computation from N = 800 to 
N = 12000. For all these runs the front stopped at the 
same position xq ~ —2, 73, with a width w ~ 1, 80. This 
phenomenon is also called "pinning of the front" 0. We 
will analyze it in detail below. 




FIG. 4: Interaction of the front with a wide "tanh" defect. 
The parameters are s; = 0, 3, s r = 1 (jump of 0, 7) and d = 10. 



When the defect is narrow, the speed and width be- 
come non monotonic. The results are shown in figure [5] 
Again the results have been checked by varying system- 
atically the spatial resolution. 
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FIG. 5: Interaction of the front with a narrow "tanh" defect. 
The parameters are s; = 0, 3, s r = 1 (jump of 0, 7) and d = 
0,1. 



As previously the dip in the velocity indicates that 
increasing the amplitude we will obtain the pinning of 
the front. This indeed happens. We will show in section 
IVII that we can infere a simple criterion for the pinning 
of the front. 



IV. COLLECTIVE COORDINATE ANALYSIS 



FIG. 3: Pinning of the front when it hits a large defect. The 
parameters are so = 0, 3, si = 7 and d = 0, 3. 



B. Tanh defect 



The numerical results in the previous section can be 
understood in a simple way. Assume that the solution 
is a function of the reduced variable z — x — ct where 
the speed c is a slowly varying function of time. Then 
reporting u into the partial differential equation (JXJ) and 
integrating we get 



We now consider a different type of defect which con- 
nects two different values s; on the left and s r on the 
right. 



. , s r — si ( x\ 
s{x) = si H — [I + tanh -j 



(8) 



We observe similar results as for the gaussian defect. For 
a "wide" defect the front moves adiabatically as shown 
in figure HI It speeds up and its width decreases as s 
increases from left to right. 



-+oo 

<■ l U'(x — ct)dx 



oo 
-t-oo 



U"(x - ct)dx + 



+ 00 



s(x)R(u(x,t))dx (9) 



Consider a front going from 1 (to the left) to 0. Since 
the front is flat away from the defect we have 



U"(z)dz = 0, / U'(z)dz = -1 



4 



so that we obtain 



+00 



s(x)R(u(x, t))dx. 



(10) 



This result shows that the speed decreases (resp. in- 
creases) when the front reaches a region where s is smaller 
(resp. larger). This is not completely correct because the 
width of the front also depends on s as shown in ((3J). 

A more general approach is then to assume that the 
front keeps its general form but that its center xo and 
width w become modulated 



u{x, t) = U 



x {t) 



w(t) 



U(z) 



(11) 



At this point we do not specify the form of U(z). This 
method is often called "variational approach" in the 
context of Hamiltonian partial differential equations Q 
which are derived from a Lagrangian density. Here how- 
ever there is no variational principle so we need to gen- 
erate two conservation laws to obtain the evolution of 
xq and w. The first one is the original partial differen- 
tial equation (p}. The second conservation law can be 
obtained by multiplying the equation by u 



uu t — uu xx + s(x)uR(u). 



(12) 



In principle we could have also multiplied by u x but this 
would give the wrong evolution if s is singular. After 
replacing the expression (|11|) into the two partial differ- 
ential equations ([T|). (IT2l . we get the evolutions of the 
front position xq and width w 



J J +w J s ( wz + x o)R = 0, 

x' J UU'+w' J UU'z+w J s{wz+x )UR = i J U' 2 , 

(13) 

where the terms fU', JU'z,fUU', JlIU'z are inte- 
grals with respect to the variable z and therefore are 
numbers. The integrals 



s{wz + Xq)R 



s(wz + xq)UR = 



s(wz + xo)R(U(z))dz, 



+00 



s(wz + xo)U(z)R(U(z))dz, 



depend on xq and w. They yield the source terms for 
the differential equations. All the details are given in 
Appendix B. These ordinary differential equations are 
very general and can be transposed to different types of 
nonlinearities. The only assumptions are that the defect 
s(x) is localized and that the front keeps its functional 
profile. 



We consider the simplest ansatz jBJ which is suitable 
for the Zeldovich equation because it is an exact solution 
when s is constant. Then ((15)) reduce to 



+00 



s(wz + xo)R(U(z))dz, 



1 

3w 



+00 



s(wz + x )(l - 2U(z))R(U(z))dz. 

Note that we still have made no assumptions on the non- 
linearity R. 

We will consider three main situations, a defect that 
varies on a scale longer than the width of the front and 
two sharp defects represented respectively as a Dirac dis- 
tribution and a Heaviside function. For a defect that 
varies on a scale longer than w, s(xq + wz) w s(xq) so 
that s goes out of the integral, modifying the equations 
subsequently 



1 - 2a 
2 

3u> 6 



ws(xo), 



1 w 
w = -ts(xo). 



(15) 



These equations justify the notion of a local speed and 
width of the front. When s(x) = s is constant, we re- 
cover the results ((4|, ((5]). Another remark is that here 
the defect can be extended. The only condition that we 
require is that the front has reached it's equilibrium state 
before reaching the defect region. 

Note that the approximation can be sharpened by use 
of a Taylor expansion : s(xq + wz) w s(xq) + wzs'(xq). 
Then the equations take the form 



. 1 w . . (1 — 2a)w 2 ,. 
W = ^-6 3{Xo) + 2 S( * 0) - 



(16) 



We now assume a sharp "bump-like" defect described 
by s(x) — a + f3S(x) where 8(x) is the Dirac distribu- 
tion centered at x = 0. The system of equations is now 
reduced to 



in = aw 



1 - 2a 



+ /3RIU 



-Xq 

w 



The source terms are given by 



R{U{z)) = -e 



(17) 



— 1 + a + ae z 



(l-2U(z))R(U(z)) = e/(l-e*) 



They are shown in Fig. |B1 



(l + e z ) 3 ' 

1 + a + ae z 



(18) 
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FIG. 6: Plot of the source terms R(U(— z)) (left panel) and 
(1 — 2U{— z)) R(U( — z)) (right panel) for a Dirac distribution 
defect as a function of the reduced variable z, for a chosen 
parameter a = 0.3. 

Note how the two sources terms R(U(—z)) and 
(1 — 2U{— z)) R(U(—z)) are not symmetric. Because 
R(U(—z)) is largely positive a defect will cause an ac- 
celeration of the kink for /3 > 0. But because there is 
also a negative part, the pinning of the front becomes 
possible both for (3 > and (3 < 0, see section PVTI for de- 
tails. Notice also how (1 — 2Z7(— zj) R(U(—z)) is largely 
negative so that the width of the front decreases as it hits 
the defect for /? > 0. 

The other sharp defect that we study is s(x) = a + 
j3H{x) where H (x) is the Heaviside function. The system 
of equations is now reduced to 



The general remarks on the acceleration and width re- 
duction of the front still apply. 

To summarize, we have obtained fairly simple ordinary 
differential equations for the evolution of the front posi- 
tion xq and width w for the Zeldovich equation. We will 
see in the next section that these equations yield very 
good approximations to the solution of the partial differ- 
ential equation. 



V. COMPARISON BETWEEN THE FULL 
MODEL AND THE REDUCED MODEL 

To establish the validity of the reduced model, it is 
important to compare its solutions to the ones of the 
partial differential equation. As discussed in the previous 
section, we classify the defects s(x) as wide or narrow 
depending whether w/d <C 1 or w/d> 1 where w is the 
initial width of the front and d is the width of the defect 
as defined in formulas ([7]) and ©. 



w 



aw 
1 

3w 



1 - 2a 



(3w 



W 



r+oo 
Pw I (1 



Here the source terms arc 



R(U(z))dz, 

- 2U{z))R{U{z))dz. 

(19) 



A. Adiabatic case 

When the defect is wide, equations (fTSj) give evolutions 
of the front position xq and width w that are close to the 
fits obtained from the solutions of the partial differential 
equation . The plots are shown in Fig. [5] for the gaussian 
defect of Fig. Q]and the tanh defect of Fig. @] 



+oo 



(1 



R(U(z))dz 
2U{z)) R{U{z))dz 



They are plotted in Fig. [71 






FIG. 8: Plots (xo, w) for the partial differential equation ([2]) in 
continuous line (red online) and the simple collective variable 
equations (115[l in dashed line (blue online). The left panel 
(resp. right panel) corresponds to a gaussian (resp. tanh) 
defect. The defects are shown at the bottom of each panel. 



FIG. 7: Plot of the source terms f+™ R(U(z))dz (left panel) 

and f_ °° (1 — 2U(z)) R(U ( z))dz (right panel) for a Heaviside 
distribution defect as a function of the reduced variable y, for 
a chosen parameter a — 0.3. 



The collective variable estimates can be improved by 
calculating numerically the integrals in equation (|13[) . 
We compare in Fig. |H] the results for the partial differ- 
ential equation (f5]) and for the solutions of (IT51) for the 
"tanh" defect. As can be seen the agreement is excellent. 
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FIG. 9: Plot (xo, w) for the partial differential equation ([2]) in 
continuous line (red online) and the collective variable equa- 
tions (fl3|) in dashed line (blue online) for the "tanh" defect. 



B. Narrow defect: pinning 

We now consider that the defect width d is smaller than 
the front width w. Then the collective variable equations 
can be reduced as shown in the previous section. We 
approximate the gaussian s(x) = sq + s\ cxp ("^"J by a 
Dirac delta function s(x) = a + /3S(x). We choose a = sq 
and j3 = si\f2ird so that 

The integrals associated to the defect are then equal. We 
tested the validity of this approximation and found that 
the solutions agree to about 1 % when w > lOd. 

For such narrow defects, the collective variable equa- 
tions P7|) are less accurate than for a wide defect. They 
do provide however the qualitative behavior, in partic- 
ular the pinning of the front. Fig. [TU] shows the plots 
for a gaussian defect drawn at the bottom. The value 
of s at infinity is sq = 0, 3. The left panel corresponds 
to a small amplitude si = 0, 6. The right panel is for 
a much larger amplitude si = 5 causing the pinning of 
the front. For this particular plot we show the speed x' Q 
computed using a centered difference for both the partial 
differential equation and the collective variables. They 
are multiplied by 10 in the plot for clarity. The defect 
has been divided by 10. 




X Xo 



FIG. 10: Plots (xo,w) for the partial differential equation 
(f2]) in continuous line (red online) and the collective variable 
equations (|17p in dashed line (blue online) for a narrow gaus- 
sian defect shown at the bottom of the plots. The left panel 
corresponds to si = 0, 6 and d = 0,3. The right panel corre- 
sponds to si = 5. There the speed x' is also reported. 



The partial differential equation and the collective vari- 
ables agree well for the width w as a function of xo . The 
speed x' is not so well approximated but it goes to zero 
for a pinning position xq. 

For narrow "tanh" defects, we reduce the collective 
variable equations to the ones for a Heaviside defect (fll?)) 
as shown above. The agreement between the curves 
(xq,w) for the partial differential equation ([2]) and the 
collective variable equations (fTU)) is good as shown in Fig. 

HH 




xo 

FIG. 11: Plot (xo, w) for the partial differential equation (f2]) in 
continuous line (red online) and the collective variable equa- 
tions (|19|l in dashed line (blue online) for a narrow "tanh" 
defect shown at the bottom of the plot. The parameters are 
si — 0, 3, s r = 1 -corresponding to a jump of 0, 7- and d = 0, 1. 

As we have seen in section 11111 the front can get pinned 
when si is large enough and this is predicted also by 
the the collective variable model. Such an example is 
shown in Fig. [12] for a large and narrow defect where 
,s r = 8, d — 0,1 shown scaled by 0.1 at the bottom of the 
plot. As previously the speed x' estimated using finite 
differences has been plotted in the graph to show pinning. 




FIG. 12: Plot (xo, w) for the partial differential equation ([2]) in 
continuous line (red online) and the collective variable equa- 
tions (|19[) in dashed line (blue online) for a narrow "tanh" 
defect shown at the bottom of the plot. The speed x' Q is also 
reported using the same color scheme. 



VI. ESTIMATES FOR FRONT PINNING AND 
DEFECT TOPOGRAPHY 

We now illustrate how the reduced model (fT^ji can be 
used. Its main advantage is that the parameters appear 
explicitely and therefore their influence can be under- 
stood. A first direct application is a simple criterion for 
the front pinning. 
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From the equation for the Dirac delta function defect, 
we can obtain a rough estimate of the strength f3 of the 
defect necessary to stop the front. The evolution of the 
front from (JTTJ) shows that the front can stop x' Q = 0, w' = 
when 





fl-2a 


aw 


V 2 


1 


w 


3w 



+ (3R \U 
P\l-2U 



w 
-xo 





R | U 



-x 



(21) 

Therefore the front can stop when the two terms on the 
right hand side of the first equation balance each other. 
In other words we need that 



1 - 2a 



s$ -/?min(i?). 



So there exist a threshold for the pinning : 

'1 - 2a s 



= — fi c min(-R). 



(22) 



The minimum of R(U(z)) as a function of z and its ar- 
gument z m i n can be computed exactly, 



= In 



1 



min(i?) = R(U(z min )) = - 



a 2 (l + r)(a + r) 
(1 + a + r) 3 ' 



(23) 
(24) 



where r = Vl — a + a 2 . Plugging j3 — j3 c into the second 
equation of (p?Tj) and assuming = ^ s - = z m i n , we get an 
expression for the width at pinning w. Using this we also 
get the front stopping position x$. For the critical (3, the 
position and width at pinning are given by 



x 



w 



•3(1 -2a) 



l + 3(l-2a)(i=*g 



(25) 



The critical (3 is obtained from (|2"2")l with the above value 
of w. It reads 



(1 - 2a) (1 + a + rf 



a*(l + r)(a + r)^l + 3(l-2a)(±=§£ 

These values are in good agreement with the values ob- 
tained numerically in the previous section, see Fig. |3] 
Precisely, with a = 0, 3, and a — 0, 3 too, we find /3 C rs 
5,88, and the parameters at pinning xq ps —3,47, w ps 
1, 89. Numerically, we found j3 c ps 6, 18, xq pa —2, 73, w ~ 
1,80. This j3 c corresponds to Si ~ 7,80 because we set 
the width of the defect d = 0, 1. Such a criterion can also 



be infered in the case of a Heaviside defect. A pinning of 
the front is also possible with j3 < (as long as j3 > —a 
so that s(x) remains positive). 

Another direct application of the collective variable 
model is that we can obtain the defect topography s(x) 
from the observation of the front position xq and width 
w. We illustrate this on the example of a wide gaus- 
sian defect of the form {7]) with ([2]) and the parameters 
so = 0, 6, si = 0, 3, d = 10. The position and width of 
the front are estimated using the least square fit on the 
solution of the partial differential equation ((2J) . From the 
collective variable equations in the adiabatic case (TTBl we 
get 



s(x ) 



1 — 2a w 



(27) 



Using a centered difference approximation for the time 
derivative, we obtain an estimate of s(xo). This estimate 
is compared to the "real" s(x) in Fig. [TBI 




FIG. 13: Defect topography s(x) estimation from the solution 
of the partial differential equation. The estimated defect to- 
pography s(x) using (|27|) is shown in dashed line (blue online) 
while the "real" s(x) is shown in continuous line (red online). 



As can be seen the agreement is very good. This ap- 
proach can then be extended to other types of defects for 
physical or biological applications. 

These two examples illustrate the power of these re- 
duced models. The role of the parameters is very easy 
to understand. One can easily solve the inverse problem 
of estimating these parameters from measurements. An- 
other extension of this could be to control the front using 
these reduced equations. 



VII. CONCLUSION 

We have analyzed numerically the interaction of Zel- 
dovich reaction-diffusion front with a localized defect. 
The stability of the front for different types of defects 
suggested that it has the form of a generalized traveling 
wave with a time dependant position and width. Us- 
ing conservation laws we obtained ordinary differential 
equations for these collective variables. We further re- 
duced these models for the cases of an adiabatic defector 
a sharp "gaussian" or "tanh" defect. For these three 



8 



cases the position and width obtained by fitting the nu- these reduced models can be used to predict the pinning 
merical solution agree very well with the solutions of the of the front on a large defect or to estimate the defect to- 
collective variable equations. Finally we illustrated how pography from a time-series of front positions and width. 
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Appendix A: Numerical procedures 

The basis of the method is to discretize the spatial 
part of the operator and keep the temporal part as such. 
We thereby transform the partial differential equation 
into a system of ordinary differential equations. This 
method allows to increase the precision of the approxi- 
mation in time and space independently and easily. We 
choose as space discretisation the finite volume approx- 
imation where the operator is integrated over reference 
volumes. The value of the function is assumed constant 
in each volume. As solver for the system of differential 
equations, we use the Runge-Kutta method of order 4- 
5 introduced by Dormand and Prince which enables to 
control the local error by varying the time-step. This 
has been implemented as the Fortran code DOPRI5 by 
Hairer and Norsett (loj . 

We introduce reference volumes Vk (interval) whose 
centers are at discretisation points Xk — x m i n + h/2 + 



(k - l)h, l^fc^n with h = (x max - x min )/n 

X k + Xk-l Xk+l + x k 



, 1 < k < n. 



For a fixed t, we assume u(x,t) to be constant on each 
volume Vk, u(xk,t) = Uk- Integrating over Vk, we obtain 
for 1 < fc < n 

u fc+ i +u k -i -2u k f .. . 

u k = jn h s k R(Uk), (Al) 

The boundary points k = 1, n are such that one satisfies 
the homogeneous Neumann boundary conditions at x — 

•Evnin , ^max • ^Ve have 



u 2 - Ml 
ui = — — h R(ui)si, 



-u„ + u n -i 
U„ = — h H{u n )s r , 



h 2 



(A2) 



(A3) 



The fit of the solutions using the front profile is done 
using a least square method. For each time t we introduce 
an error function 



1 " 

E(x , w) = - y"Vtti(t) - u k (x l , t)Y 



(A4) 



where the ui are the values of u calculated numerically 
and 



u k (x,t) 



1 



exp 



( x-x (t) \ 

V Wi ) 



is the exact kink solution. The expression E is mini- 
mized usin g th e Polak-Ribiere combination of line min- 
imisations The fit of the solutions is done on the 
region such that 0.01 < u k < 0.99. The initial guesses 
for xq and w are estimated from u w 0.25 and u w 0.75. 
The method converges in about 20 iterations and the 
value of the energy is small, min E < 10~ 4 . Therefore 
the fit is good. 



Appendix B: Derivation of the collective variable 
equations 



The traveling wave is 
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so that we have the relations for the partial derivatives 
1 / , w'{x — xq)\ , fx — Xq' 



Ui 



1 , fx - x 
u x = —U 



w 



U' 



w 



The second conservation law (|T2j) will yield the evolu- 
tion of w. Proceeding as above we get 



w \ w 
1 tth fx- x o 



w z \ w 



, (B2) 
(B3) 
(B4) 



The first equation is obtained by integrating (JT|) with 
respect to x 

/+oo r + oo 

u t dx = [uxjx^Loo + / s(x)R(u)dx 
-oo J — oo 

The term [ux^^L^ is zero because the front is flat away 
from the defect. Introducing the partial derivatives we 
get 



U' ( ° }dx I ^ (X - Xq) fx- Xq 



/+oo r+oo 
U(z)U'(z)dz + w' U(z)U'(z)zdz 
-oo J —oo 

/ + OO 1 f+OG 

s(wz+x )U(z)R(U(z))dz = - / U' 2 (z)dz 



(B7) 



w \ w 

f + OO 



s(x)R 



dx (B5) 



Note that we only assumed U = 1 at — oo and U = 
at oo and that [/' is even. These assumptions are very 
\general. In particular we have made no restrictions on 
ihfecreaction term R. 

For the Zeldovich reaction term R(u) — u(l — u)(u — a) 
it is natural to assume that 



We compute the integrals by making the change of vari- 
ables X = W Z + Xq 

r+co r+co 

-x' Q U'(z)dz-w' / U'(z)zdz 



-Hoc 



U(z) 



1 



— w s(wz + xo)R(U(z))dz 



1 + exp z 

Then the integrals not involving s can be computed and 
we obtain the final result. 



Assuming that the front goes from 1 to and that U'(z) 
is even we get the final result 



x' / U'(z)dz + w' / U'(z)zdz 



+ w s(wz + x o )R(U(z))dz = (B6) 



x' Q = w / s(wz + xo)R(U(z))dz 



1 r +oc 

w' = — + w / s(wz + x )(l - 2U(z))R(U(z))dz 



(B8) 



